setwd("code/")
source("insurance_source_score_1_prep.R")

setwd("code/")
source("insurance_source_score_3_functions.R")

lmers_market_acafavor <- plot_slopes_return_lmer_list(
    outcome="FAVOR >= 3",
    propensity = "MARKET",
  the_subset=subset(
    sdat.under65.2,
    ## we were originally going to use different income cutoffs for ACAST 1 vs 0
    (
      (ACAST==1 & INCOME>40)
      | (ACAST==0 & INCOME>40)
    )
  ),
  plot_title="Exchange Score\nrespondents under 65, income over 40K",
  pdf_name="figureA8_pscore-by-survey-20190208-income-cutoffs-40k-with-weights-lm.pdf",
  no_pdf = FALSE
)

## median: 50
with(subset(
    sdat.under65.2,
    INCOME>40 & DATN > 60 & MARKET
  ), median(AGEN))



lmers_market_partyid <- plot_slopes_return_lmer_list(
    outcome="DEM",
    the_subset=subset(
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Exchange Score\nrespondents under 65, income over 40k",
    pdf_name="figureA10_pscore-by-survey-20190104-partyid-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="blue", no_pdf=FALSE
)




mod_for_pid_diff <- lm(
    FAVOR >= 3 ~ GOP,
    data=subset(
        sdat.under65.2,
         INCOME>40 & PERIOD > 8.5 & (IND==0)
    ),
    weights=WEIGHT## weights are survey weights
)
summary(mod_for_pid_diff)
coef(mod_for_pid_diff)[2] # used for uninsured comparison in text of paper


#### Figure 5
pdf(
    "../figs/figure5_selfinsure-uninsured-pscore-by-survey-20190104-income-cutoffs-40k-with-weights-lm.pdf",
    height=3.15, width=10
)
par(mar=c(6, 4, 4, 3), mfrow=c(1,2))
lmers_selfinsured_acafavor <- plot_slopes_return_lmer_list(
    propensity="SELFINSURE==1",
    outcome="FAVOR>=3",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Self-insured Score\nrespondents under 65, income over 40k",
    pdf_name="figure5_selfinsure-pscore-by-survey-20190104-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="ACA favorability",
    .the_prefix=the_prefix, no_pdf=T,
    color="orange4"
)

    par(xpd = F)
lmers_uninsured_acafavor <- plot_slopes_return_lmer_list(
    propensity="COVERED==0",
    outcome="FAVOR>=3",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Uninsured Score\nrespondents under 65, income over 40k",
    pdf_name="figure5_uninsured-pscore-by-survey-20190104-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="ACA favorability",
    .the_prefix=the_prefix, no_pdf=T,
    color=gray(0.2)## "purple"
)

dev.off()

lmers_uninsured_acafavor <- plot_slopes_return_lmer_list(
    propensity="COVERED==0",
    outcome="FAVOR>=3",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Uninsured Score\nrespondents under 65, income over 40k",
    pdf_name="uninsured-pscore-by-survey-20190104-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="ACA favorability",
    .the_prefix=the_prefix, no_pdf=F,
    color=gray(0.2)## "purple"
)

coef(summary(lmers_uninsured_acafavor[["lmer_out_all"]]))[4] / coef(mod_for_pid_diff)[2]



lmers_uninsured_partyid <- plot_slopes_return_lmer_list(
    propensity="COVERED==0",
    outcome="DEM",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Uninsured Score\nrespondents under 65, income over 40k",
    pdf_name="uninsured-pscore-by-survey-20190104-partyid-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="blue"
)





lmers_selfinsured_acafavor <- plot_slopes_return_lmer_list(
    propensity="SELFINSURE==1",
    outcome="FAVOR>=3",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Self-insured Score\nrespondents under 65, income over 40k",
    pdf_name="selfinsure-pscore-by-survey-20190104-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="ACA favorability",
    .the_prefix=the_prefix, no_pdf=F,
    color="orange4"
)



stargazer(
    lmers_market_acafavor[["propensity_out_lm"]],
    lmers_uninsured_acafavor[["propensity_out_lm"]],
    lmers_selfinsured_acafavor[["propensity_out_lm"]],
    omit=c("as.factor\\(ST\\)"),
    omit.stat=c("ser","adj.rsq", "rsq"),
    table.layout = "ts", float=F, digits=2,
    out=paste0(
        the_prefix, "/figs/",
        "tableA17_insurance_propensity_lm_models.tex"
    )
)


exponentiate <- function(x) exp(x)

## lm_simple is actually a poisson model
stargazer(
    lmers_market_acafavor[["propensity_out_lm_simple"]],
    lmers_uninsured_acafavor[["propensity_out_lm_simple"]],
    lmers_selfinsured_acafavor[["propensity_out_lm_simple"]],
    ci=TRUE,
    apply.coef=exponentiate, ## apply.ci=exponentiate,
    omit.stat=c("ser","adj.rsq", "rsq"),
    table.layout = "ts", float=F, digits=2,
    out=paste0(
        the_prefix, "/figs/",
        "tableA21_insurance_source_score_interpretable_coefficients.tex"
    )
)

stargazer(
    lmers_market_acafavor[[1]], lmers_uninsured_acafavor[[1]], lmers_selfinsured_acafavor[[1]],
    covariate.labels=c(
        "Propensity Score", "Post-Implementation",
        "Post-Implementation x Propensity Score",
        "Intercept"
    ),
    ## omit="Constant",
    omit.stat=c("ser","f","adj.rsq", "rsq"),
    table.layout = "ts", float=F, digits=2,
    out=paste0(
        the_prefix, "/figs/",
        "tableA20_insurance_source_score_multilevel_models.tex"
    )
)


print(
    xtable(
        sapply(
            data.frame(
                `Market propensity`=lmers_market_acafavor[["out_data"]]$PSCORE3_orig,
                `Uninsured propensity`=lmers_uninsured_acafavor[["out_data"]]$PSCORE3_orig,
                `Self-insured propensity`=lmers_selfinsured_acafavor[["out_data"]]$PSCORE3_orig
            ),
            function(x)
                data.frame(
                    Mean=mean(x, na.rm=T),
                    SD=sd(x, na.rm=T),
                    Min=min(x, na.rm=T),
                    Max=max(x, na.rm=T)
                )
        )
    ),
    digits=2,
    file=paste0(
        the_prefix, "/figs/",
        "tableA18_insurance_source_score_summary_statistics.tex"
    )
)


## well outside bounds without scaling
lmers_selfinsured_partyid <- plot_slopes_return_lmer_list(
    propensity="SELFINSURE==1",
    outcome="DEM",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Self-insured Score\nrespondents under 65, income over 40k",
    pdf_name="selfinsure-pscore-by-survey-20190104-partyid-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="orange3",
    color2="blue"
)


## summary(lmers_selfinsured_partyid[[1]])

lmers_market_partyid_control_acafavor <- plot_slopes_return_lmer_list(
    outcome="DEM",
    the_subset=subset(
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Exchange Score\nrespondents under 65, income over 40k",
    pdf_name="pscore-by-survey-20190104-partyid-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="blue",
    control_acafavor=TRUE, no_pdf=TRUE
)

lmers_uninsured_partyid_control_acafavor <- plot_slopes_return_lmer_list(
    propensity="COVERED==0",
    outcome="DEM",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Uninsured Score\nrespondents under 65, income over 40k",
    pdf_name="uninsured-pscore-by-survey-20190104-partyid-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="blue",
    control_acafavor=TRUE, no_pdf=TRUE
)

lmers_selfinsured_partyid_control_acafavor <- plot_slopes_return_lmer_list(
    propensity="SELFINSURE==1",
    outcome="DEM",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Self-insured Score\nrespondents under 65, income over 40k",
    pdf_name="selfinsure-pscore-by-survey-20190104-partyid-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="orange3",
    color2="blue",
    control_acafavor=TRUE, no_pdf=TRUE
)


lmers_market_independent <- plot_slopes_return_lmer_list(
    outcome="PID==2",
    the_subset=subset(
        sdat.under65.2,
        ## EMPLINSURE==0
        ## &
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Exchange Score\nrespondents under 65, income over 40k",
    pdf_name="figureA11_pscore-by-survey-20190104-independent-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Independent",
    .the_prefix=the_prefix,
    color="green4",
    ## color2="green4",
    no_pdf=FALSE
)

lmers_uninsured_independent <- plot_slopes_return_lmer_list(
    propensity="COVERED==0",
    outcome="PID==2",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Uninsured Score\nrespondents under 65, income over 40k",
    pdf_name="uninsured-pscore-by-survey-20190104-independent-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Independent",
    .the_prefix=the_prefix,
    color="green4",
    no_pdf=TRUE
)

lmers_selfinsured_independent <- plot_slopes_return_lmer_list(
    propensity="SELFINSURE==1",
    outcome="PID==2",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Self-insured Score\nrespondents under 65, income over 40k",
    pdf_name="selfinsure-pscore-by-survey-20190104-independent-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Independent",
    .the_prefix=the_prefix,
    color="orange3",
    color2="green4",
    no_pdf=TRUE
)


lmers_market_independent_control_acafavor <- plot_slopes_return_lmer_list(
    outcome="PID==2",
    the_subset=subset(
        sdat.under65.2,
        ## EMPLINSURE==0
        ## &
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Exchange Score\nrespondents under 65, income over 40k",
    pdf_name="pscore-by-survey-20190104-independent-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="blue",
    control_acafavor=TRUE, no_pdf=TRUE
)

lmers_uninsured_independent_control_acafavor <- plot_slopes_return_lmer_list(
    propensity="COVERED==0",
    outcome="PID==2",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Uninsured Score\nrespondents under 65, income over 40k",
    pdf_name="uninsured-pscore-by-survey-20190104-independent-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="blue",
    control_acafavor=TRUE, no_pdf=TRUE
)

lmers_selfinsured_independent_control_acafavor <- plot_slopes_return_lmer_list(
    propensity="SELFINSURE==1",
    outcome="PID==2",
    the_subset=subset(                  # on training set only
        sdat.under65.2,
        (
            (ACAST==1 & INCOME>40)
            | (ACAST==0 & INCOME>40)
        )
    ),
    plot_title="Self-insured Score\nrespondents under 65, income over 40k",
    pdf_name="selfinsure-pscore-by-survey-20190104-independent-control-acafavor-income-cutoffs-40k-with-weights-lm.pdf",
    y_axis="Democrat",
    .the_prefix=the_prefix,
    color="orange3",
    color2="blue",
    control_acafavor=TRUE, no_pdf=TRUE
)


cormatrix_pred <- cor(
    as.matrix(
        cbind(
            lmers_market_acafavor[["out_data"]]$PSCORE3,
            lmers_uninsured_acafavor[["out_data"]]$PSCORE3,
            lmers_selfinsured_acafavor[["out_data"]]$PSCORE3
        )
    ), use="complete"                   #RETIRED missing for 6k
)


print(
    xtable(
        cbind(cormatrix_pred),
        digits=2
    ), include.rownames=F,
    file=paste0(
        the_prefix, "/figs/",
        "tableA19_insurance_source_score_correlation_matrix.tex"
    )
)







## must be order self-insured (1) uninsured (2)

plot_comparisons(
    mod_obj1=lmers_selfinsured_partyid,
    mod_obj2=lmers_uninsured_partyid,
    pdf_name="figureA10_selfinsure-vs-uninsure-pscore-by-survey-20190104-partyid-with-weights-lm.pdf",
    y_axis="Democrat",
    color="blue3",
    color2=gray(0.4)
)

plot_comparisons(
    mod_obj1=lmers_selfinsured_independent,
    mod_obj2=lmers_uninsured_independent,
    pdf_name="figureA11_selfinsure-vs-uninsure-pscore-by-survey-20190104-independent-with-weights-lm.pdf",
    y_axis="Independent",
    color="green4",
    color2=gray(0.4)
)


model_list <- list(
    lmers_market_partyid[["lmer_out_all"]],
    lmers_market_partyid_control_acafavor[["lmer_out_all"]],
    lmers_market_independent[["lmer_out_all"]],
    lmers_market_independent_control_acafavor[["lmer_out_all"]],
    lmers_uninsured_partyid[["lmer_out_all"]],
    lmers_uninsured_partyid_control_acafavor[["lmer_out_all"]],
    lmers_uninsured_independent[["lmer_out_all"]],
    lmers_uninsured_independent_control_acafavor[["lmer_out_all"]],
    lmers_selfinsured_partyid[["lmer_out_all"]],
    lmers_selfinsured_partyid_control_acafavor[["lmer_out_all"]],
    lmers_selfinsured_independent[["lmer_out_all"]],
    lmers_selfinsured_independent_control_acafavor[["lmer_out_all"]]
)
##
pdf(
    paste0(
        the_prefix, "/figs/",
        "figureA9_all-partyid-pscore-by-survey-20190224-with-weights-lm.pdf"
    ),
    width=6, height=3.5
    )
par(mar=c(4,5,1,1))
plot(
    x=0,
    y=0,
    xlim=c(1,12), ylim=c(-1, 1.2),
    bty="n",
    type="n",
    xaxt="n",
    ylab="Change in Partisanship\n2014 on vs. 2013 and earlier",
    xlab="Insurance Source Scores",
    cex.lab=1
)
abline(
    h=0, lty=2
)
abline(
    v=c(4.5, 8.5), lty=3
)
j <- as.vector(rbind(
    seq(1.25, 11.25, 2),
    seq(1.6, 11.6, 2)
    ))
for (i in 1:length(model_list)) {
    points(
        x=j[i],
        y=summary(model_list[[i]])$coef["PSCORE3:I(PERIOD > 8.5)TRUE","Estimate"],
        col=ifelse(i %in% c(3:4, 7:8, 11:12), "green4", "blue3"),
        pch=ifelse(i %in% seq(1, 11, 2), 16, 15),
        cex=ifelse(i %in% seq(1, 11, 2), 1.6, 0.9)
    )
    segments(
        x0=j[i],
        x1=j[i],
        y0=summary(model_list[[i]])$coef["PSCORE3:I(PERIOD > 8.5)TRUE","Estimate"] -
                                   1.96 * summary(model_list[[i]])$coef["PSCORE3:I(PERIOD > 8.5)TRUE","Std. Error"],
        y1=summary(model_list[[i]])$coef["PSCORE3:I(PERIOD > 8.5)TRUE","Estimate"] +
                                   1.96 * summary(model_list[[i]])$coef["PSCORE3:I(PERIOD > 8.5)TRUE","Std. Error"],
        col=ifelse(i %in% c(3:4, 7:8, 11:12), "green4", "blue3"),
        lwd=ifelse(i %in% seq(1, 11, 2), 4, 1.5)
    )
}
axis(
    side=1,
    at=seq(2.5, 10.5, 4),
    labels=c("Exchanges", "Uninsured", "Self-insured"),
    tick=F,
    cex.axis=1.2
)
legend(
    x=3.9,
    y=1.42,
    legend=c("DV: Shift in Party ID\n(D=1; I=0.5; R=0)","DV: Shift to Independence"),
    pch=c(16, 16),
    col=c("blue3","green4"),
    bty="n",## title="DV:",
    pt.cex=1.3
)
legend(
    x=9,
    y=0.8,
    legend=c("Control ACA\nfavorability"),
    pch=c(15),
    col=c("black"),
    bty="n"
)
dev.off()
